if (runTime.outputTime())
  {
    volScalarField epsilonEq
      (
       IOobject
       (
	"epsilonEq",
	runTime.timeName(),
	mesh,
	IOobject::NO_READ,
	IOobject::AUTO_WRITE
	),
       sqrt((2.0/3.0)*magSqr(dev(epsilon)))
       );

    Info<< "Max epsilonEq = " << max(epsilonEq).value()
	<< endl;

    volScalarField sigmaEq
      (
       IOobject
       (
	"sigmaEq",
	runTime.timeName(),
	mesh,
	IOobject::NO_READ,
	IOobject::AUTO_WRITE
	),
       sqrt((3.0/2.0)*magSqr(dev(sigma)))
       );
    
    Info<< "Max sigmaEq = " << max(sigmaEq).value()
	<< endl;
    
    //- Calculate Cauchy stress
    volTensorField F = I + gradU;
    volScalarField J = det(F);

    //- update density
    rho = rho/J;

    volSymmTensorField sigmaCauchy
      (
       IOobject
       (
	"sigmaCauchy",
	runTime.timeName(),
	mesh,
	IOobject::NO_READ,
	IOobject::AUTO_WRITE
	),
       (1/J) * symm(F.T() & sigma & F)
       );

    //- Cauchy von Mises stress
    volScalarField sigmaCauchyEq
      (
       IOobject
       (
	"sigmaCauchyEq",
	runTime.timeName(),
	mesh,
	IOobject::NO_READ,
	IOobject::AUTO_WRITE
	),
       sqrt((3.0/2.0)*magSqr(dev(sigmaCauchy)))
       );
    
    Info<< "Max sigmaCauchyEq = " << max(sigmaCauchyEq).value()
	<< endl;

    //- write boundary forces
    //- integrate (sigma2PK & F) over reference area
    //- which is equivalent to integrating sigmaCauchy
    //- over the deformed area
    Info << nl;
    forAll(mesh.boundary(), patchi)
      {
	Info << "Patch " << mesh.boundary()[patchi].name() << endl;
	tensorField F = I + gradU.boundaryField()[patchi];
	vectorField totalForce = mesh.Sf().boundaryField()[patchi] & (sigma.boundaryField()[patchi] & F);

	vector force = sum( totalForce );
	Info << "\ttotal force is " << force << " N" << endl;
	
	tensorField Finv = inv(F);
	vectorField nCurrent = Finv & n.boundaryField()[patchi];
	nCurrent /= mag(nCurrent);
	scalar normalForce = sum( nCurrent & totalForce );
	Info << "\tnormal force is " << normalForce << " N" << endl;
	scalar shearForce = mag(sum( (I - sqr(nCurrent)) & totalForce ));
	Info << "\tshear force is " << shearForce << " N" << endl << endl;;
      }


    //- move mesh for visualisation and move it back after writing
//     vectorField oldPoints = mesh.allPoints();
//     #include "moveMeshLeastSquares.H"
    
    runTime.write();

    //- move mesh back
//     mesh.movePoints(oldPoints);
  }
